setwd("~/Desktop/Appendix")
pkgs <- c( "systemfit", "dummies", "sandwich", "lmtest", "stargazer","countrycode", "assertthat", "readr", "fBasics", "texreg")
invisible(lapply(pkgs, function(x) if(!is.element(x, installed.packages()[, 1])) install.packages(x, repos = c(CRAN = "http://cran.rstudio.com"))))
invisible(lapply(pkgs, require, character.only = TRUE))
options(stringsAsFactors = FALSE)
options(scipen=999)

dat <- read.csv("~/Desktop/Appendix/Data Analysis/final_dat.csv")

#####Summary Statistics #### 
### This can take over 24 hours to complete ####### 
temp <- dat[,c("l.access_points", "workerr.mean2", "latentmean",  "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.fdiflows_UNCTAD", "l.cescr_ratify",  "l.cat_ratify", "l.conflict", "l.under_BCG", "l.polity")]

sum_stat_table <- round(t(basicStats(temp)[c("Mean", "Stdev", "Median", "Minimum", "Maximum", "NAs"),]), 4)
rownames(sum_stat_table) <- c("Access Points", "Workers' Rights", "Physical Integrity Rights", "log(GDP per capita)", "log(Population)", "Youth Bulge", "Trade", "Oil", "FDI", "CESCR Ratify", "CAT Ratify", "Conflict", "IMF Conditionality", "Polity")

###### Multiple Imputation ##### 
library(Amelia)
library(systemfit)

dat <- dat[, c(1,2,12,17, 28:31,33:39, 41:47,50,51)]
dat <- dat[dat$year > 1979, ]



set.seed(082592)
impute_dat <- amelia(x=dat, 
                     m=5, 
                     ts="year", 
                     cs="ccode", 
                     polytime=3, 
                     intercs=TRUE, 
                     empri= 0.01 * nrow(dat)) 


#----------------------------------------------------------------------------------------------------------------------------------------------------------------#

for(i in 1:5){  
  
  
  
  library(dummies)
  country.dummies <- as.data.frame(dummy(impute_dat$imputations[[i]]$ccode))
  
  
  
  colnames(country.dummies) <- c("ccode20" ,"ccode31" ,"ccode41" ,"ccode42" ,"ccode51" ,"ccode52" ,"ccode53" ,"ccode54" ,"ccode55" ,"ccode56" ,"ccode57" ,"ccode58" ,"ccode60" ,"ccode70" ,"ccode80" ,"ccode90" ,"ccode91" ,"ccode92" ,"ccode93" ,"ccode94" ,"ccode95" ,"ccode100" ,"ccode101" ,"ccode110" ,"ccode115" ,"ccode130" ,"ccode135" ,"ccode140" ,"ccode145" ,"ccode155" ,"ccode160" ,"ccode165" ,"ccode200","ccode205" ,"ccode210" ,"ccode211" ,"ccode212" ,"ccode220" ,"ccode223" ,"ccode225" ,"ccode230" ,"ccode232" ,"ccode235" ,"ccode255" ,"ccode290" ,"ccode305" ,"ccode310" ,"ccode315" ,"ccode316" ,"ccode317" ,"ccode325" ,"ccode331" ,"ccode338" ,"ccode339" ,"ccode343" ,"ccode344" ,"ccode349" ,"ccode350" ,"ccode352" ,"ccode355" ,"ccode359" ,"ccode360" ,"ccode365" ,"ccode366" ,"ccode367" ,"ccode368" ,"ccode369" ,"ccode371" ,"ccode375" ,"ccode380" ,"ccode385" ,"ccode390" ,"ccode395" ,"ccode402" ,"ccode403" ,"ccode404" ,"ccode432" ,"ccode433" ,"ccode434" ,"ccode436" ,"ccode437" ,"ccode451" ,"ccode452" ,"ccode475" ,"ccode482" ,"ccode484" ,"ccode500" ,"ccode501" ,"ccode516" ,"ccode551" ,"ccode553" ,"ccode560" ,"ccode565" ,"ccode570" ,"ccode580" ,"ccode581" ,"ccode590" ,"ccode625" ,"ccode640" ,"ccode666" ,"ccode712" ,"ccode713" ,"ccode732" ,"ccode740" ,"ccode750" ,"ccode770" ,"ccode771" ,"ccode780" ,"ccode790" ,"ccode800" ,"ccode840" ,"ccode850" ,"ccode900" ,"ccode910" ,"ccode920" ,"ccode935" ,"ccode940" ,"ccode946" ,"ccode970" ,"ccode983" ,"ccode986" ,"ccode987")  
  
  impute_dat$imputations[[i]] <- cbind(impute_dat$imputations[[i]], country.dummies)
  
} 


latentmean.model2 <- latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade + l.log.oil  + l.conflict + l.cat_ratify + l.polity 

workerr.model2 <-  workerr.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade  + l.fdi.permillion + l.conflict  +  l.imf_con + l.cescr_ratify + l.polity



latentmean.model2.fe <- latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade + l.log.oil  + l.conflict + l.cat_ratify + l.polity + ccode41  +  ccode42  +  ccode51  +  ccode80 + ccode90 + ccode91 + ccode92 + ccode93 + ccode94 + ccode95 + ccode100 + ccode101 + ccode110 +  ccode130 + ccode140 + ccode145  + ccode160  + ccode339 +  ccode355 + ccode360 +  ccode402 +  ccode432  + ccode434 + ccode436 +  ccode451  + ccode475 + ccode482 + ccode551 + ccode553 + ccode560 + ccode581 + ccode590 +  ccode640 +  ccode712 +  ccode750 + ccode771 + ccode780 + ccode790 + ccode800 + ccode840 


workerr.model2.fe <-  workerr.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade  + l.fdi.permillion + l.conflict  +  l.imf_con + l.cescr_ratify + l.polity +  ccode41  +  ccode42  +  ccode51  +  ccode80 + ccode90 + ccode91 + ccode92 + ccode93 + ccode94 + ccode95 + ccode100 + ccode101 + ccode110 +  ccode130 + ccode140 + ccode145  + ccode160  + ccode339 +  ccode355 + ccode360 +  ccode402 +  ccode432  + ccode434 + ccode436 +  ccode451  + ccode475 + ccode482 + ccode551 + ccode553 + ccode560 + ccode581 + ccode590 +  ccode640 +  ccode712 +  ccode750 + ccode771 + ccode780 + ccode790 + ccode800 + ccode840 




#----------------------------------------------------------------------------------------------------------------------------------------------------------------# 

impute_sur <- function(model1, model2, impute_dat){ 
  
  
  b.out<- NULL 
  se.out<- NULL      
  
  for(i in 1:impute_dat$m){ 
    
    #Access Point Variable  
    fitsur <- systemfit(list(surfit.latentmean2=model1, surfit.workerr2=model2), method= "SUR", data=impute_dat$imputations[[i]])
    
    
    b.out<- rbind(b.out, summary(fitsur)$coefficients[,1])  
    se.out <- rbind(se.out, summary(fitsur)$coefficients[,2])
    
  }    
  
  
  combine.results <- mi.meld(q=b.out, se=se.out)
  
  results<- matrix(NA, length(combine.results$q.mi), 4)      
  results[,1]<- combine.results$q.mi
  results[,2]<- combine.results$se.mi 
  results[,3]<- results[,1]/results[,2] 
  results[,4]<- 2*pnorm(-abs(results[,3]))
  rownames(results) <- colnames(b.out)
  
  round(results,8) 
  
  
  
}

impute_sur(latentmean.model2, workerr.model2, impute_dat)
impute_sur(latentmean.model2.fe, workerr.model2.fe, impute_dat)


#########Robustness Check - Fixed Effects###########




library(bucky) 
papermod_pir0 <- mi.eval(lm(latentmean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity, data=impute_dat), robust=TRUE)
summary(papermod_pir0)
papermod_pir1 <- mi.eval(lm(latentmean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity +  factor(ccode), data=impute_dat), robust=TRUE)
summary(papermod_pir1)
papermod_pir2 <-mi.eval(lm(latentmean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity +  factor(year), data=impute_dat), robust=TRUE)
summary(papermod_pir2)
papermod_pir3 <- mi.eval(lm(latentmean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity +  factor(ccode) + factor(year), data=impute_dat), robust=TRUE)
summary(papermod_pir3)



papermod_wr0 <- mi.eval(lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=impute_dat), robust=TRUE)
summary(papermod_wr0)
papermod_wr1 <- mi.eval(lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil + factor(ccode), data=impute_dat), robust=TRUE)
summary(papermod_wr1)
papermod_wr2 <- mi.eval(lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil + factor(year), data=impute_dat), robust=TRUE)
summary(papermod_wr2)
papermod_wr3 <- mi.eval(lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil + factor(ccode) + factor(year), data=impute_dat), robust=TRUE)
summary(papermod_wr3)



######## Conditional Analysis ######## 

dat<- final_dat[complete.cases(final_dat[,c("latentmean", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil","l.conflict", "l.cat_ratify", "l.polity", "l.fdiflows_UNCTAD", "l.lji_LS", "l.cescr_ratify", "l.under_BCG", "workerr.mean2", "l.CSOConsult")]), ]

library(ggplot2)

ggplot(dat, aes(x=l.CSOConsult)) + 
  geom_density() +
  theme_light() + 
  labs(x="CSO Consult")


m2.r<- lm(formula = latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity + l.CSOConsult + l.access_points*l.CSOConsult, data = dat)

m3.r<- lm(formula = latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity + l.CSOConsult + l.access_points*l.CSOConsult + factor(ccode), data = dat)


robust.se <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- model$rank
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  rcse.se <- coeftest(model, rcse.cov)
  return(rcse.se)
}

se2.r <- robust.se(m2.r, cluster=dat$ccode)

se3.r <- robust.se(m3.r, cluster=dat$ccode)

library(interplot)
interplot(m= m2.r, var1 = "l.access_points", var2="l.CSOConsult", hist=TRUE)  + 
  xlab("CSO Consult") + 
  ylab("Estimated Coefficient for \nAccess Points") + 
  theme_light()+ 
  geom_hline(yintercept = 0, linetype="dashed")

interplot(m= m2.r, var1 = "l.CSOConsult", var2="l.access_points", hist=TRUE)  + 
  ylab("Estimated Coefficient for \n CSO Consult ") + 
  xlab("Access Points") + 
  theme_light()+ 
  geom_hline(yintercept = 0, linetype="dashed")


###### Control for Judicial Independence #### 

papermod_pir <- mi.eval(lm(latentmean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity + l.lji_LS, data=impute_dat), robust=TRUE)
papermod_wr <- mi.eval(lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil + l.lji_LS, data=impute_dat), robust=TRUE)                      
 

######## Alternative Dependent Variables ######### 

hr_scores <- read.csv("~/Desktop/Appendix/Data Analysis/Source Data/HumanRightsProtectionScores_v2.04.csv")
colnames(hr_scores) <- c("year", "ciri.name", "ccode", "disap", "kill", "polpris", "tort", "amnesty", "pts", "hathaway", "itt", "genocide", "rummel", "massive_repression", "executions","killings", "ciri", "latentmean", "latentsd") 

workerr_sources <-read.csv("~/Desktop/Appendix/latent_workerr/Sources/workerr_sources.csv")

appdat<- merge(dat, hr_scores, by=c("ccode", "year"), all.y=FALSE)
appdat<- merge(appdat, workerr_sources, by=c("ccode", "year"), all.y=FALSE)
 

robust.se <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- model$rank
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  rcse.se <- coeftest(model, rcse.cov)
  return(rcse.se)
}


cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "pts")]), ]
mod6 <- lm(pts  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity, data=cmp.appdat) 
se6 <- robust.se(mod6, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "ciri")]), ]
mod14 <- lm(ciri  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity  , data=cmp.appdat) 
se14 <- robust.se(mod14, cmp.appdat$ccode)[,2] 


cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod16 <- lm( workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se16 <- robust.se(mod16, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "minimum.wage", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod17 <-lm( minimum.wage  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se17 <- robust.se(mod17, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "association.govt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod18 <- lm(association.govt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se18 <- robust.se(mod18, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "association.mrkt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod19 <- lm(association.mrkt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se19 <- robust.se(mod19, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "colbargain.govt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod20 <- lm(colbargain.govt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se20<- robust.se(mod20, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "colbargain.mrkt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod21 <- lm(colbargain.mrkt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se21<- robust.se(mod21, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "strike.govt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod22 <- lm(strike.govt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se22<- robust.se(mod22, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "strike.mrkt", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod23 <- lm(strike.mrkt  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se23<- robust.se(mod23, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "ciri_worker", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod24 <- lm(ciri_worker  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se24<- robust.se(mod24, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "LaborRightsPos", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod25 <- lm(LaborRightsPos  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se25<- robust.se(mod25, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "LawPos", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod26 <- lm(LawPos  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se26 <- robust.se(mod26, cmp.appdat$ccode)[,2] 

cmp.appdat <- appdat[complete.cases(appdat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cescr_ratify", "l.polity", "PracticePos", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
mod27 <- lm(PracticePos  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
se27 <- robust.se(mod27, cmp.appdat$ccode)[,2] 

###### Worker Rights Latent Variable Country Trends ####### 

library(ggplot2)
library(countrycode)
#breaks it up into multiple sections to not overload memory #### 
temp<- final_dat[final_dat$ccode < 101, ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p1<- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)

print(p1) 

temp<- final_dat[final_dat$ccode %in% c(101:309), ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p2<- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)

print(p2)


temp<- final_dat[final_dat$ccode %in% c(310:380), ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p3 <- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)


print(p3) 

temp<- final_dat[final_dat$ccode %in% c(385:565), ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p4<- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)

print(p4)

temp<- final_dat[final_dat$ccode %in% c(570:910), ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p5 <- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)

print(p5)

temp<- final_dat[final_dat$ccode > 910, ]
temp$name <- countrycode(temp$ccode, "cown", "country.name")

p6 <- ggplot(temp, aes(x=year, y=workerr.mean)) + 
  geom_line() + 
  facet_wrap(~name, ncol=3)

print(p6)


temp_work <- workerr_latent[workerr_latent$year==2002, ]
temp_work1 <- temp_work[temp_work$ccode %in% c(2,20,31,40,41,42,51,52,53,54,55,56,57,58,60,70,80,90,91,92,93,94,95,100,101,110,115,130,135,140,145,155,160,165,200,205,210,211,212,220,223,225,230,232,235,255,290,305,310,315,316,317,325),]
temp_work2 <- temp_work[temp_work$ccode %in% c(331,338,339,343,344,349,350,352,355,359,360,365,366,367,368,369,371,375,380,385,390,395,402,404,432,433,434,436,451,452,475,482,484,500,501,516,520,551,553,560,565,570,580,581,590,625,640,660,666,712,713,732,740,750,770,771,775,780,790,800,812,840,850,900,910,920,935,940,946,970,983,986,9870),]
p7 <- ggplot(temp_work1, aes(x=country, y=workerr.mean)) + 
  geom_point() + 
  coord_flip() +
  geom_errorbar(aes(ymin=workerr.mean - workerr.sd, ymax=workerr.mean+workerr.sd))

print(p7)

p8 <- ggplot(temp_work2, aes(x=country, y=workerr.mean)) + 
  geom_point() + 
  coord_flip() +
  geom_errorbar(aes(ymin=workerr.mean - workerr.sd, ymax=workerr.mean+workerr.sd))

print(p8)

############ Uncertainty Models ########### 


hr_scores <- read.csv("~/Desktop/Appendix/HumanRightsProtectionScores_v2.04.csv")
colnames(hr_scores) <- c("year", "ciri.name", "ccode", "disap", "kill", "polpris", "tort", "amnesty", "pts", "hathaway", "itt", "genocide", "rummel", "massive_repression", "executions","killings", "ciri", "latentmean", "latentsd") 
hr_scores <- hr_scores[, c("year", "ccode", "latentsd")] 
assert_that(!anyDuplicated(hr_scores[, c("ccode", "year")]))


dat<- merge(final_dat, hr_scores, by=c("ccode","year"), all.x=TRUE, all.y=FALSE) 
colnames(dat)

dat <- dat[complete.cases(dat[,c("latentmean", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil","l.conflict", "l.cat_ratify", "l.polity", "l.fdiflows_UNCTAD", "l.lji_LS", "l.cescr_ratify", "l.under_BCG", "workerr.mean2")]), ]


dat$id <-as.numeric(as.factor(dat$ccode))




library(R2jags)

hrlatent <- dat$latentmean
hrlatent_var <- dat$latentsd
workerr.mean2 <- dat$workerr.mean2
workerr.mean2_var <- dat$workerr.sd2
l.access_points <- dat$l.access_points
l.log.realgdppccurrent <- dat$l.log.realgdppccurrent
l.log.pop <- dat$l.log.pop
l.log.oil <- dat$l.log.oil
l.ythblgap <- dat$l.ythblgap
l.trade_WDI <- dat$l.trade_WDI
l.conflict <- dat$l.conflict
l.cat_ratify <- dat$l.cat_ratify 
l.polity <- dat$l.polity
l.fdiflows_UNCTAD <- dat$l.fdiflows_UNCTAD
l.under_BCG <- dat$l.under_BCG 
l.cescr_ratify <- dat$l.cescr_ratify
l.CSOConsult <- dat$l.CSOConsult
id <- dat$id
year <- dat$year 

N <- length(dat$ccode) 

bayes.dat <- list("hrlatent", "hrlatent_var", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.conflict", "l.cat_ratify", "l.log.oil", 
                  "l.polity", "N")

bayes.params <- c("alpha", "b1", "b2", "b3", "b4", "b5", "b6", "b7", "b8", "b9", "tau")

bayes.inits1 <- list("alpha"= -.5, "b1"= -.5, "b2"= -.5, "b3"= -.5, "b4"= -.5, "b5" = -.5, "b6"=-.5, "b7"= -.5, "b8"= -.5, "b9"=-.5)
bayes.inits2 <- list("alpha"= .5, "b1"= .5, "b2"= .5, "b3"= .5, "b4"= .5, "b5" = .5, "b6"= .5, "b7"= .5, "b8"= .5, "b9" = .5)

bayes.inits <- list(bayes.inits1, bayes.inits2)

cat('data {
  for (i in 1:N) {
    
    hr_var[i] <- (1/hrlatent_var[i])  
    
    hrlatent_est[i]~dnorm(hrlatent[i], hr_var[i]) 
    
    #work_var[i] <- (1/workerr.mean2_var[i]) 
    
    #work_est[i]~dnorm(workerr.mean2[i], work_var[i]) 
      
  }
}
model{
  for(i in 1:N){
    mu[i] <-  alpha
    + b1*l.access_points[i]
    + b2*l.log.realgdppccurrent[i] 
    + b3*l.log.pop[i]
    + b4*l.ythblgap[i]
    + b5*l.trade_WDI[i]
    + b6*l.conflict[i]
    + b7*l.cat_ratify[i]
    + b8*l.polity[i]
    + b9*l.log.oil[i]
    hrlatent_est[i]~dnorm(mu[i], tau)
  }
  
  
    alpha~dnorm(0,.01)
  
  
  b1~dnorm(0,.01)
  b2~dnorm(0,.01)
  b3~dnorm(0,.01)
  b4~dnorm(0,.01)
  b5~dnorm(0,.01)
  b6~dnorm(0,.01)
  b7~dnorm(0,.01)
  b8~dnorm(0,.01)
  b9~dnorm(0,.01)

  tau~dgamma(.1, 1)
}', file={PIRlatent_bayes <- tempfile()})


pirfit.1 <- jags(data=bayes.dat, inits=bayes.inits,
                 parameters=bayes.params, n.chains=2, n.iter=100000, n.burnin=20000,
                 model.file=PIRlatent_bayes)

pirfit.1.upd <- update(pirfit.1, n.iter=300000, n.burnin=20000)






bayes.dat2 <- list("id", "hrlatent", "hrlatent_var", "l.access_points",
                   "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.conflict", "l.cat_ratify", "l.polity", "l.log.oil", "l.CSOConsult", "N")

bayes.params2 <- c("alpha", "b1", "b2", "b3", "b4", "b5", "b6", "b7", "b8", "b9", "b10", "b11", "tau")

bayes.inits1.2 <- list("alpha"= -.5, "b1"= -.5, "b2"= -.5, "b3"= -.5, "b4"= -.5, "b5" = -.5, "b6"=-.5, "b7"= -.5, "b8"= -.5, "b9"=-.5, "b10"= -.5, "b11" = -.5)
bayes.inits2.2 <- list("alpha"= .5, "b1"= .5, "b2"= .5, "b3"= .5, "b4"= .5, "b5" = .5, "b6"= .5, "b7"= .5, "b8"= .5, "b9" = .5, "b10" = .5, "b11" = .5) 

bayes.inits2 <- list(bayes.inits1.2, bayes.inits2.2)

cat('data {
  for (i in 1:N) {
    
    hr_var[i] <- (1/hrlatent_var[i])  
    
    hrlatent_est[i]~dnorm(hrlatent[i], hr_var[i]) 
    
    #work_var[i] <- (1/workerr.mean2_var[i]) 
    
    #work_est[i]~dnorm(workerr.mean2[i], work_var[i]) 
      
  }
}
model{
  for(i in 1:N){
    mu[i] <-  alpha
    + b1*l.access_points[i]
    + b2*l.log.realgdppccurrent[i] 
    + b3*l.log.pop[i]
    + b4*l.ythblgap[i]
    + b5*l.trade_WDI[i]
    + b6*l.conflict[i]
    + b7*l.cat_ratify[i]
    + b8*l.polity[i]
    + b9*l.log.oil[i]
    + b10*l.CSOConsult[i]
    + b11*(l.CSOConsult[i] * l.access_points[i])
    hrlatent_est[i]~dnorm(mu[i], tau)
  }
  
  
    alpha~dnorm(0,.01)
  
  
  b1~dnorm(0,.01)
  b2~dnorm(0,.01)
  b3~dnorm(0,.01)
  b4~dnorm(0,.01)
  b5~dnorm(0,.01)
  b6~dnorm(0,.01)
  b7~dnorm(0,.01)
  b8~dnorm(0,.01)
  b9~dnorm(0,.01)
  b10~dnorm(0,.01)
  b11~dnorm(0,.01)

  tau~dgamma(.1, 1)
}', file={PIRlatent_bayes2 <- tempfile()})

pirfit.2 <- jags(data=bayes.dat2, inits=bayes.inits2,
                 parameters=bayes.params2, n.chains=2, n.iter=100000, n.burnin=20000,
                 model.file=PIRlatent_bayes2)

pirfit.2.upd <- update(pirfit.2, n.iter=300000, n.burnin=20000)




bayes.dat3 <- list("id", "workerr.mean2", "workerr.mean2_var", "l.access_points",
                   "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.conflict", 
                   "l.polity", "l.fdiflows_UNCTAD", "l.under_BCG", "l.cescr_ratify", "l.log.oil", "N")

bayes.params3 <- c("alpha", "b1", "b2", "b3", "b4", "b5", "b6", "b12", "b8", "b9", "b13", "b14", "tau")

bayes.inits1.3 <- list("alpha"= -.5, "b1"= -.5, "b2"= -.5, "b3"= -.5, "b4"= -.5, "b5" = -.5, "b6"=-.5, "b12"= -.5, "b8"= -.5, "b9"=-.5, "b13"= -.5, "b14" = -.5)
bayes.inits2.3 <- list("alpha"= .5, "b1"= .5, "b2"= .5, "b3"= .5, "b4"= .5, "b5" = .5, "b6"= .5, "b12"= .5, "b8"= .5, "b9" = .5, "b13" = .5, "b14" = .5) 

bayes.inits3 <- list(bayes.inits1.3, bayes.inits2.3)

cat('data {
  for (i in 1:N) {
    
    #hr_var[i] <- (1/hrlatent_var[i])  
    
    #hrlatent_est[i]~dnorm(hrlatent[i], hr_var[i]) 
    
    work_var[i] <- (1/workerr.mean2_var[i]) 
    
    work_est[i]~dnorm(workerr.mean2[i], work_var[i]) 
      
  }
}
model{
  for(i in 1:N){
    mu[i] <-  alpha
    + b1*l.access_points[i]
    + b2*l.log.realgdppccurrent[i] 
    + b3*l.log.pop[i]
    + b4*l.ythblgap[i]
    + b5*l.trade_WDI[i]
    + b6*l.conflict[i]
    + b12*l.cescr_ratify[i]
    + b8*l.polity[i]
    + b9*l.log.oil[i]
    + b13*l.fdiflows_UNCTAD[i]
    + b14*l.under_BCG[i]
    work_est[i]~dnorm(mu[i], tau)
  }
  
  
    alpha~dnorm(0,.01)
  
  
  b1~dnorm(0,.01)
  b2~dnorm(0,.01)
  b3~dnorm(0,.01)
  b4~dnorm(0,.01)
  b5~dnorm(0,.01)
  b6~dnorm(0,.01)
  b12~dnorm(0,.01)
  b8~dnorm(0,.01)
  b9~dnorm(0,.01)
  b13~dnorm(0,.01)
  b14~dnorm(0,.01)

  tau~dgamma(.1, 1)
}', file={PIRlatent_bayes3 <- tempfile()})

pirfit.3 <- jags(data=bayes.dat3, inits=bayes.inits3,
                 parameters=bayes.params3, n.chains=2, n.iter=100000, n.burnin=20000,
                 model.file=PIRlatent_bayes3)

pirfit.3.upd <- update(pirfit.3, n.iter=300000, n.burnin=20000)










library(texreg)
library(BayesPostEst)

mcmcReg(list(pirfit.1.upd, pirfit.2.upd, pirfit.3.upd),
        custom.coef.map = list('alpha' = "(Intercept)", 
                               'b1' = "Access Points", 
                               'b2' = "log(GDP per capita)",
                               'b3' = "log(Population)",
                               'b4' = "Youth Bulge",
                               'b5' = "Trade",
                               'b6' = "Conflict",
                               'b7' = "CAT Ratify",
                               'b8' = "Polity",
                               'b9' = "Oil", 
                               'b10'= "CSO Consult", 
                               'b11'= "Access Points x CSO Consult",
                               'b12' = "CESCR Ratify", 
                               'b13' = "FDI", 
                               'b14' = "IMF Conditionality"), 
        custom.model.names = c("PIR", "PIR with interaction", "Worker's Rights"),
        ci=.9,
        digits=3) 


library(texreg)
library(BayesPostEst)

mcmcReg(list(pirfit.1.upd, pirfit.2.upd, pirfit.3.upd),
        custom.coef.map = list('alpha' = "(Intercept)", 
                               'b1' = "Access Points", 
                               'b2' = "log(GDP per capita)",
                               'b3' = "log(Population)",
                               'b4' = "Youth Bulge",
                               'b5' = "Trade",
                               'b6' = "Conflict",
                               'b7' = "CAT Ratify",
                               'b8' = "Polity",
                               'b9' = "Oil", 
                               'b10'= "CSO Consult", 
                               'b11'= "Access Points x CSO Consult",
                               'b12' = "CESCR Ratify", 
                               'b13' = "FDI", 
                               'b14' = "IMF Conditionality"), 
        custom.model.names = c("PIR", "PIR with interaction", "Worker's Rights"),
        ci=.8,
        digits=3) 



######### Alternative Specifications of the Latent Worker Variable #########

library(R2jags)
dat <-read.csv("~/Desktop/Appendix/latent_workerr/Sources/workerr_sources.csv")

#remove data where everything is unknown 
datsplit <- split(dat, dat$ccode)
dat <- NULL
for (i in 1:length(datsplit)) {
  yvars <- datsplit[[i]][, c("minimum.wage", "association.govt", "association.mrkt", "colbargain.govt", "colbargain.mrkt", "strike.govt", "strike.mrkt", "ciri_worker", "LaborRightsPos", "LawPos", "PracticePos", "uri_l", "uri_p", "sri_l", "sri_p", "wri_l", "wri_p")] 
  first <- as.numeric(which(cumsum(!apply(is.na(yvars), 1, all))==1))
  last <- nrow(datsplit[[i]])-as.numeric(which(cumsum(rev(!apply(is.na(yvars), 1, all)))==1))+1
  if (length(first) > 1) {first <- min(first)}
  if (length(last) > 1) {last <- min(last)}
  if (length(first) > 0) {
    dat <- rbind(dat, datsplit[[i]][first:last, ])
  }
}


#Define Vectors of the data matrix for JAGS
minimum.wage <- dat$minimum.wage
association.govt <- dat$association.govt 
association.mrkt <- dat$association.mrkt 
colbargain.govt <- dat$colbargain.govt 
colbargain.mrkt <- dat$colbargain.mrkt 
strike.govt <- dat$strike.govt 
strike.mrkt <- dat$strike.mrkt 
mosley.law <- dat$LawPos
workr.union.law <-dat$uri_l
workr.substantive.law <- dat$sri_l
N <- nrow(dat) 

#Read in the Civil Liberties data for JAGS
workerr.data <- list("minimum.wage", "association.govt", "association.mrkt", "colbargain.govt", "colbargain.mrkt", "strike.govt", "strike.mrkt", "mosley.law", "workr.union.law",  "workr.substantive.law", "N")

#Define the parameters interested in:
workerr.params <- c("workerr", "b", "tau")


#'cat' the model because the path isn't being found
cat('model {
for (i in 1:N) {
  minimum.wage[i]~dnorm(mu1[i], tau[1])
  mu1[i]<-b[1]*workerr[i]
  
  association.mrkt[i]~dnorm(mu2[i], tau[2])
  mu2[i]<-b[2]*workerr[i]
  
  association.govt[i]~dnorm(mu3[i], tau[3])
  mu3[i]<-b[3]*workerr[i]
  
  strike.mrkt[i]~dnorm(mu4[i], tau[4])
  mu4[i]<-b[4]*workerr[i]

  strike.govt[i]~dnorm(mu5[i], tau[5])
  mu5[i]<-b[5]*workerr[i]

  colbargain.mrkt[i] ~ dnorm(mu6[i], tau[6])
  mu6[i]<-b[6]*workerr[i] 

  colbargain.govt[i] ~ dnorm(mu7[i], tau[7])
  mu7[i]<-b[7]*workerr[i] 

  mosley.law[i] ~ dnorm(mu8[i], tau[8])
  mu8[i]<-b[8]*workerr[i]

  workr.union.law[i] ~ dnorm(mu9[i], tau[9]) 
  mu9[i] <- b[9]*workerr[i]

  workr.substantive.law[i] ~ dnorm(mu10[i], tau[10]) 
  mu10[i] <- b[10]*workerr[i]
  

}

for (i in 1:N) {
  workerr[i]~dnorm(0, 1)
}

for (i in 1:10) {
  b[i]~dnorm(0, .1)T(0,)
}

for (i in 1:10) {
  tau[i]~dgamma(1, 1)
}
}', file={workerr <- tempfile()})


#Fit the model in JAGS
workerrfit <- jags(data=workerr.data, inits=NULL, workerr.params, n.chains = 2, n.iter = 10000, n.burnin=1000, model.file = workerr)


#calculate the mean for each column
workerr.law.mean <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, mean)
workerr.law.sd <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, sd)

#create df
workerr.law.lats <- cbind(workerr.law.mean, workerr.law.sd)
dat <- cbind(dat, workerr.law.lats)








ciri_worker <- dat$ciri_worker 
mosley.practice <- dat$PracticePos
workr.union.practice <-dat$uri_p
workr.substantive.practice <-dat$sri_p
N <- nrow(dat) 

#Read in the Civil Liberties data for JAGS
workerr.data <- list("ciri_worker", "mosley.practice", "workr.union.practice",  "workr.substantive.practice", "N")

#Define the parameters interested in:
workerr.params <- c("workerr", "b", "tau")


#'cat' the model because the path isn't being found
cat('model {
for (i in 1:N) {

  mosley.practice[i] ~ dnorm(mu1[i], tau[1])
  mu1[i]<-b[1]*workerr[i]

  ciri_worker[i] ~ dnorm(mu2[i], tau[2]) 
  mu2[i] <- b[2]*workerr[i]
  
  
  workr.union.practice[i] ~ dnorm(mu3[i], tau[3]) 
  mu3[i] <- b[3]*workerr[i]
 
  workr.substantive.practice[i] ~ dnorm(mu4[i], tau[4]) 
  mu4[i] <- b[4]*workerr[i]

}

for (i in 1:N) {
  workerr[i]~dnorm(0, 1)
}

for (i in 1:4) {
  b[i]~dnorm(0, .1)T(0,)
}

for (i in 1:4) {
  tau[i]~dgamma(1, 1)
}
}', file={workerr <- tempfile()})


#Fit the model in JAGS
workerrfit <- jags(data=workerr.data, inits=NULL, workerr.params, n.chains = 2, n.iter = 10000, n.burnin=1000, model.file = workerr)


#calculate the mean for each column
workerr.practice.mean <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, mean)
workerr.practice.sd <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, sd)

#create df
workerr.practice.lats <- cbind(workerr.practice.mean, workerr.practice.sd)
dat <- cbind(dat, workerr.practice.lats)



final_dat <- read.csv("~/Desktop/Appendix/Data Analysis/final_dat.csv")


test_dat <- merge(final_dat, dat, by=c("ccode", "year"), all.x=TRUE, all.y=FALSE)



cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod1 <- lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse1<- robust.se(testmod1, cmp.appdat$ccode)[,2] 


cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod2 <- lm(workerr.law.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse2<- robust.se(testmod2, cmp.appdat$ccode)[,2] 


cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod3 <- lm(workerr.practice.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse3<- robust.se(testmod3, cmp.appdat$ccode)[,2] 


########## Alternative Specifications of the Latent Worker Variable ########### 

library(R2jags)
dat <-read.csv("~/Desktop/Appendix/latent_workerr/Sources/workerr_sources.csv")

#remove data where everything is unknown 
datsplit <- split(dat, dat$ccode)
dat <- NULL
for (i in 1:length(datsplit)) {
  yvars <- datsplit[[i]][, c("minimum.wage", "association.govt", "association.mrkt", "colbargain.govt", "colbargain.mrkt", "strike.govt", "strike.mrkt", "ciri_worker", "LaborRightsPos", "LawPos", "PracticePos", "uri_l", "uri_p", "sri_l", "sri_p", "wri_l", "wri_p")] 
  first <- as.numeric(which(cumsum(!apply(is.na(yvars), 1, all))==1))
  last <- nrow(datsplit[[i]])-as.numeric(which(cumsum(rev(!apply(is.na(yvars), 1, all)))==1))+1
  if (length(first) > 1) {first <- min(first)}
  if (length(last) > 1) {last <- min(last)}
  if (length(first) > 0) {
    dat <- rbind(dat, datsplit[[i]][first:last, ])
  }
}


#Define Vectors of the data matrix for JAGS

association.govt <- dat$association.govt 
association.mrkt <- dat$association.mrkt 
colbargain.govt <- dat$colbargain.govt 
colbargain.mrkt <- dat$colbargain.mrkt 
strike.govt <- dat$strike.govt 
strike.mrkt <- dat$strike.mrkt 
mosley.law <- dat$LawPos
workr.union.law <-dat$uri_l
mosley.practice <- dat$PracticePos
workr.union.practice <-dat$uri_p
mosley <- dat$LaborRightsPos


minimum.wage <- dat$minimum.wage
workr.substantive.law <- dat$sri_l
workr.substantive.practice <-dat$sri_p

N <- nrow(dat) 

#Read in the Civil Liberties data for JAGS
workerr.data <- list("association.govt", "association.mrkt", "colbargain.govt", "colbargain.mrkt", "strike.govt", "strike.mrkt", "mosley.law", "mosley.practice", "mosley", "workr.union.law", "workr.union.practice", "N")

#Define the parameters interested in:
workerr.params <- c("workerr", "b", "tau")


#'cat' the model because the path isn't being found
cat('model {
for (i in 1:N) {
  mosley[i]~dnorm(mu1[i], tau[1])
  mu1[i]<-b[1]*workerr[i]
  
  association.mrkt[i]~dnorm(mu2[i], tau[2])
  mu2[i]<-b[2]*workerr[i]
  
  association.govt[i]~dnorm(mu3[i], tau[3])
  mu3[i]<-b[3]*workerr[i]
  
  strike.mrkt[i]~dnorm(mu4[i], tau[4])
  mu4[i]<-b[4]*workerr[i]

  strike.govt[i]~dnorm(mu5[i], tau[5])
  mu5[i]<-b[5]*workerr[i]

  colbargain.mrkt[i] ~ dnorm(mu6[i], tau[6])
  mu6[i]<-b[6]*workerr[i] 

  colbargain.govt[i] ~ dnorm(mu7[i], tau[7])
  mu7[i]<-b[7]*workerr[i] 

  mosley.law[i] ~ dnorm(mu8[i], tau[8])
  mu8[i]<-b[8]*workerr[i]

  workr.union.law[i] ~ dnorm(mu9[i], tau[9]) 
  mu9[i] <- b[9]*workerr[i]

  workr.union.practice[i] ~ dnorm(mu10[i], tau[10]) 
  mu10[i] <- b[10]*workerr[i]
  
  mosley.practice[i] ~ dnorm(mu11[i], tau[11]) 
  mu11[i] <- b[11]*workerr[i]

}

for (i in 1:N) {
  workerr[i]~dnorm(0, 1)
}

for (i in 1:11) {
  b[i]~dnorm(0, .1)T(0,)
}

for (i in 1:11) {
  tau[i]~dgamma(1, 1)
}
}', file={workerr <- tempfile()})


#Fit the model in JAGS
workerrfit <- jags(data=workerr.data, inits=NULL, workerr.params, n.chains = 2, n.iter = 10000, n.burnin=1000, model.file = workerr)


#calculate the mean for each column
workerr.union.mean <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, mean)
workerr.union.sd <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, sd)

#create df
workerr.union.lats <- cbind(workerr.union.mean, workerr.union.sd)
dat <- cbind(dat, workerr.union.lats)










#Read in the Civil Liberties data for JAGS
workerr.data <- list("minimum.wage", "workr.substantive.practice", "workr.substantive.law", "N")

#Define the parameters interested in:
workerr.params <- c("workerr", "b", "tau")


#'cat' the model because the path isn't being found
cat('model {
for (i in 1:N) {

  minimum.wage[i] ~ dnorm(mu1[i], tau[1])
  mu1[i]<-b[1]*workerr[i]

  workr.substantive.law[i] ~ dnorm(mu2[i], tau[2]) 
  mu2[i] <- b[2]*workerr[i]
  
  
  workr.substantive.practice[i] ~ dnorm(mu3[i], tau[3]) 
  mu3[i] <- b[3]*workerr[i]

}

for (i in 1:N) {
  workerr[i]~dnorm(0, 1)
}

for (i in 1:3) {
  b[i]~dnorm(0, .1)T(0,)
}

for (i in 1:3) {
  tau[i]~dgamma(1, 1)
}
}', file={workerr <- tempfile()})


#Fit the model in JAGS
workerrfit <- jags(data=workerr.data, inits=NULL, workerr.params, n.chains = 2, n.iter = 10000, n.burnin=1000, model.file = workerr)


#calculate the mean for each column
workerr.substantive.mean <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, mean)
workerr.substantive.sd <- apply((workerrfit$BUGSoutput$sims.list$workerr), 2, sd)

#create df
workerr.substantive.lats <- cbind(workerr.substantive.mean, workerr.substantive.sd)
dat <- cbind(dat, workerr.substantive.lats)



final_dat <- read.csv("~/Desktop/Appendix/Data Analysis/final_dat.csv")


test_dat <- merge(final_dat, dat, by=c("ccode", "year"), all.x=TRUE, all.y=FALSE)



cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod1 <- lm(workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse1<- robust.se(testmod1, cmp.appdat$ccode)[,2] 


cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod2 <- lm(workerr.union.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse2<- robust.se(testmod2, cmp.appdat$ccode)[,2] 


cmp.appdat <- test_dat[complete.cases(test_dat[, c("ccode", "year", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil", "l.conflict", "l.cat_ratify", "l.polity", "workerr.mean2", "l.fdiflows_UNCTAD", "l.under_BCG")]), ]
testmod3 <- lm(workerr.substantive.mean  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil, data=cmp.appdat) 
testse3<- robust.se(testmod3, cmp.appdat$ccode)[,2] 

